subroutine READ_ERAmet(nstep, nx, ny, nxs2i, nxs2e,        &
   & nys2i, nys2e, uwind, vwind,q,potemp,psuf, pl, dt)
integer nstep, nx, ny, nxs2i, nxs2e, nys2i, nys2e
real(8) uwind(nxs2i:nxs2e,nys2i:nys2e), vwind(nxs2i:nxs2e,nys2i:nys2e), &
&  q(nxs2i:nxs2e,nys2i:nys2e), potemp(nxs2i:nxs2e,nys2i:nys2e), &
& psuf(nxs2i:nxs2e,nys2i:nys2e), pl(nxs2i:nxs2e,nys2i:nys2e)
integer ki, i, j
real(8) dt
real temp(nx, ny,6)
real(8) rh
real(8) temp2(nxs2i:nxs2e,nys2i:nys2e,6)
character(4), dimension(6), parameter :: inp_file=(/'u10','v10','d2m','t2m','sp','tp'/)
real(8) TTD2RH,RH2Q
do ki=1,6
OPEN(37,FILE='binary_ERA/'//trim(inp_file(ki)),STATUS='UNKNOWN',FORM='UNFORMATTED', ACCESS='DIRECT',RECL=nx*ny)
read(37,rec=nstep)((temp(i,j,ki),i=1,nx),j=1,ny)
close(37)
end do
temp2(nxs2i:0,nys2i:nys2e,1:6)=999.
temp2(nx+1:nxs2e, nys2i:nys2e,1:6)=999.
temp2(1:nx, nys2i:0,1:6)=999.
temp2(1:nx, ny+1:nys2e,1:6)=999.
temp2(1:nx, 1:ny,1:6)=real(temp(1:nx,1:ny,1:6),8)
do j=1,ny
do i=1,nx
uwind(i,j)=temp2(i,j,1)
vwind(i,j)=temp2(i,j,2)
rh=TTD2RH(temp2(i,j,4),temp2(i,j,3))
q(i,j)=RH2Q(temp2(i,j,3),temp2(i,j,5)/100.,rh)
psuf(i,j)=temp2(i,j,5)
potemp(i,j)=temp2(i,j,4)*(100000./temp2(i,j,5))**0.286
pl(i,j)=dmax1(0.0, temp2(i,j,6))
!write(*,*)uwind(i,j),vwind(i,j), q(i,j),psuf(i,j),potemp(i,j),
!write(*,*) pl(i,j)
end do
end do

end subroutine READ_ERAmet

subroutine READ_ERArad(nstep, nx, ny, nxs2i, nxs2e,        &
   & nys2i, nys2e, Srad_surf,Lrad_surf,SraddirPAR_surf,SraddirIR_surf,SraddifPAR_surf,SraddifIR_surf,dt)
integer nstep, nx, ny, nxs2i, nxs2e, nys2i, nys2e
real(8) Srad_surf(nxs2i:nxs2e,nys2i:nys2e),Lrad_surf(nxs2i:nxs2e,nys2i:nys2e), &
& SraddirPAR_surf(nxs2i:nxs2e,nys2i:nys2e),SraddirIR_surf(nxs2i:nxs2e,nys2i:nys2e), &
& SraddifPAR_surf(nxs2i:nxs2e,nys2i:nys2e),SraddifIR_surf(nxs2i:nxs2e,nys2i:nys2e)
character(4), dimension(4), parameter :: inp_file=(/'tcc','ssrd','strd','tisr'/)
integer ki, i, j
real temp(nx, ny,4)
real(8) temp2(nxs2i:nxs2e,nys2i:nys2e,4)
real(8) dir, diff
real(8) dt
real(8) frac_diff
do ki=1,4
OPEN(37,FILE='binary_ERA/'//trim(inp_file(ki)),STATUS='UNKNOWN',FORM='UNFORMATTED', ACCESS='DIRECT',RECL=nx*ny)
read(37,rec=nstep)((temp(i,j,ki),i=1,nx),j=1,ny)
close(37)
end do
temp2(nxs2i:0,nys2i:nys2e,:)=999.
temp2(nx+1:nxs2e, nys2i:nys2e,:)=999.
temp2(1:nx, nys2i:0,:)=999.
temp2(1:nx, ny+1:nys2e,:)=999.
temp2(1:nx, 1:ny,1:4)=real(temp(1:nx,1:ny,1:4),8)
do j=1, ny
do i=1,nx
Srad_surf(i,j)=temp2(i,j,2)
Lrad_surf(i,j)=temp2(i,j,3)
if(temp2(i,j,4).ne.0.0)then
diff=frac_diff(temp2(i,j,2)/temp2(i,j,4))*Srad_surf(i,j)
else
diff=0.0
end if
dir=Srad_surf(i,j)-diff
SraddirPAR_surf(i,j)=0.5*dir
SraddirIR_surf(i,j)=0.5*dir
SraddifPAR_surf(i,j)=0.5*diff
SraddifIR_surf(i,j)=0.5*diff
end do
end do
end subroutine READ_ERArad
  !-----------------------------------------------------------
  !> calculate relative humidity from temperature and dew point temperature
  !-----------------------------------------------------------
  real(8) function TTD2RH(t,td)
    implicit none

    real(8), intent(in) :: t  !< temperature (in K)
    real(8), intent(in) :: td !< dew point temperature (in K)
    real(8) :: par_pres
    real(8) :: e,es

    es = par_pres(t)

    if ( t < td ) then
      e = es
   ! else if ( t > td + 25._8 ) then
    !  e = par_pres(t - 25._8)
    else
      e = par_pres(td)
    endif

    TTD2RH = max(0.03_8, min(1.23_8, e/es))
    return

  end function TTD2RH

  !-----------------------------------------------------------
  !> Saturation pressure at t temperature (K)
  !> algorithm: http://www.faqs.org/faqs/meteorology/temp-dewpoint/
  !-----------------------------------------------------------
  real(8) function par_pres(t)
    implicit none

    real(8), intent(in) :: t !< Temperature in K

    real(8) :: ct

    ct = t - 273.15_8

    if (t > 273.15_8) then
      par_pres = 6.1078_8*10**(7.5_8*ct/(ct + 237.3_8))
    else
      par_pres = 6.1078_8*10**(9.5_8*ct/(ct + 265.5_8))
    endif

  end function par_pres

  !-----------------------------------------------------------
  !> specific humidity from t, p & rh
  !-----------------------------------------------------------
  real(8) function RH2Q(t, p, rh)
    implicit none

    real(8), intent(in) :: t  !< Temperature in K
    real(8), intent(in) :: p  !< Pressure in hPa
    real(8), intent(in) :: rh !< Relative humidity in percent

    real(8) :: es, e
    real(8) :: par_pres
    real(8), parameter :: k_vap_dry = 0.62197_8

    es = par_pres(t) !min(0.5_8, par_pres(t)/p)

    e = rh * es

    RH2Q = k_vap_dry * e / (p + e*(k_vap_dry - 1._8)) !rh * a / (1._8 - rh*a*(R_vap/R_dry - 1))

!    write(*,*) 'calculating Q: es: ', es, ', p : ', p, ', rh: ', rh, ', t: ', t, ', q: ', RH2Q


  end function RH2Q
  
   real(8) function frac_diff(ss0)
   real(8) ss0
   if(ss0.le.0.22)then
   frac_diff=1-0.09*ss0
   elseif(ss0.gt.0.22.and.ss0.le.0.8)then
   frac_diff=0.9511-0.1604*ss0+4.388*ss0**2-16.638*ss0**3+12.336*ss0**4
   elseif(ss0.gt.0.8) then
   frac_diff=0.165
   end if
   end function frac_diff
